---

Homework 10 : Conditioning and stability of linear least squares¶

---

The least squares problem

$$\begin{equation} A\mathbf x = \mathbf b \end{equation}$$

where $A \in \mathcal R^{m \times n}$, $m \ge n$ has four associated "conditioning" problems, described in the table in Theorem 18.1 of TB (page 131). These are

1. Sensitivity of $\mathbf y = A\mathbf x$ to right hand side vector $\mathbf b$,

2. Sensitivity of the solution $\mathbf x$ to right hand side vector $\mathbf b$,

3. Sensitivity of $\mathbf y = A\mathbf x$ to the coefficient matrix $A$, and

4. Sensitivity of the solution $\mathbf x$ to the coefficient matrix $A$.

Problem 1¶

---

Sensitivity of $\mathbf y$ to a perturbation in $\mathbf b$.

In TB Lecture 12, the relative condition number is defined as

$$\begin{equation} \kappa = \sup_{\delta x}\left(\frac{\Vert \delta f \Vert}{\Vert f(x) \Vert}\bigg/ \frac{\Vert \delta x \Vert}{\Vert x \Vert}\right) \end{equation}$$

Problem 1(a)¶

Arguing directly from this definition, establish the condition number of $\mathbf y$ with respect to perburbations in $\mathbf b$ given by TB Lecture 18

$$\begin{equation} \kappa = \frac{1}{\cos \theta} \end{equation}$$

Hint: The input "$x$" in this problem is $\mathbf b$ and the output (or model) "$f$" is $\mathbf y$. Show geometrically that the supremum is attained with $P\delta \mathbf b = \delta \mathbf b$.

Solution¶

Since the relative condition number is defined as ; $$\begin{equation} \kappa = \sup_{\delta \mathbf x}\left(\frac{\Vert \delta \mathbf f \Vert}{\Vert \mathbf f(\mathbf x) \Vert}\bigg/ \frac{\Vert \delta \mathbf x \Vert}{\Vert \mathbf x \Vert}\right), \end{equation}$$

then the relative condition number of $\mathbf y$ with respect to perburbations in $\mathbf b$ is given by;

$$\begin{equation} \kappa = \sup_{\delta \mathbf b}\left(\frac{\Vert \delta \mathbf y \Vert}{\Vert \mathbf y \Vert}\bigg/ \frac{\Vert \delta \mathbf b \Vert}{\Vert \mathbf b \Vert}\right) \end{equation}$$

where $\mathbf y = A\mathbf x = P \mathbf b$ is the the projetion of $\mathbf b$ onto the $range(A)$, $\delta \mathbf b$ is a pertubation in $\mathbf b$ with a correponding pertubation $\delta \mathbf y$ in $\mathbf y$ (which is the projetion of $\delta \mathbf b$ onto the $range(A)$) . All this is illustrated geometrically in the image below;

From the above image we notice that $\cos (\theta) = \frac{\Vert y \Vert}{\Vert b \Vert}$ where $\theta$ is the angle between $\mathbf b$ and the $range(A)$.

The relative condition number of $\mathbf y$ with respect to perburbations in $\mathbf b$ then becomes;

$$\begin{align*} \kappa &= \sup_{\delta \mathbf b}\left(\frac{1}{\cos (\theta)}\times \frac{\Vert \delta \mathbf y \Vert}{\Vert \delta \mathbf b \Vert}\right)\\\\ &= \sup_{\delta \mathbf b}\left(\frac{1}{\cos (\theta)}\times \frac{\Vert P\delta \mathbf b \Vert}{\Vert \delta \mathbf b \Vert}\right)\\\\ &= \frac{1}{\cos (\theta)}\times \sup_{\delta \mathbf b}\left(\frac{\Vert P\delta \mathbf b \Vert}{\Vert \delta \mathbf b \Vert}\right) \end{align*}$$

But to find the $\sup_{\delta \mathbf b}\left(\frac{\Vert P\delta \mathbf b \Vert}{\Vert \delta \mathbf b \Vert}\right)$, we must maximize the ratio $\frac{\Vert P\delta \mathbf b \Vert}{\Vert \delta \mathbf b \Vert}$.

From the diagram above we see that; for $0<\beta< \pi/2$, $\Vert \delta \mathbf b \Vert$ will always be greater than $\Vert \delta \mathbf y \Vert$ and therefore to maximize $\frac{\Vert P\delta \mathbf b \Vert}{\Vert \delta \mathbf b \Vert}$, $\beta$ must be equal to zero and $\Vert \delta \mathbf b \Vert = \Vert P\delta \mathbf b \Vert$.

This means that we obtain $\sup_{\delta \mathbf b}\left(\frac{\Vert P\delta \mathbf b \Vert}{\Vert \delta \mathbf b \Vert}\right)$ when $\Vert \delta \mathbf b \Vert = \Vert P\delta \mathbf b \Vert$ and therefore the condition number of $\mathbf y$ with respect to perburbations in $\mathbf b$ becomes; $$\begin{equation} \kappa = \frac{1}{\cos \theta} \end{equation}$$

Problem 1(b)¶

For $\theta = \pi/2$, the condition number is $\infty$. Illustrate what this means by considering the least squares problem

$$\begin{equation} \begin{bmatrix} 2 \\ 1 \end{bmatrix} \begin{bmatrix} x \end{bmatrix} = \begin{bmatrix} -1 \\ 2 \end{bmatrix} \end{equation}$$

Use the results in TB 11.11 and 11.12 (page 82) to determine the projection operator $P$ for this problem. Then compute $\mathbf y = P\mathbf b$ and show that $P\mathbf b = 0$. Find a perturbation $\delta \mathbf b$ so that $P\delta \mathbf b = \delta \mathbf b = \delta \mathbf y \ne 0$. Explain what a condition number $\kappa=\infty$ might mean here. Illustrate your argument graphically.

Solution¶

From the above system we have that;

$A = \begin{bmatrix} 2 \\ 1 \end{bmatrix}$ and $\mathbf b = \begin{bmatrix} -1 \\ 2 \end{bmatrix}$

To find the the projection operator for this problem, we use $P = AA^\dagger$, but $$\begin{align*} A^\dagger &= (A^TA)A^T\\ & = \bigg( \begin{bmatrix} 2 & 1 \end{bmatrix} \begin{bmatrix} 2 \\ 1 \end{bmatrix}\bigg)^{-1} \begin{bmatrix} 2 & 1 \end{bmatrix}\\ &= \frac{1}{5}\begin{bmatrix} 2 & 1 \end{bmatrix}\\ \Rightarrow P & = \frac{1}{5}\begin{bmatrix} 2 \\ 1 \end{bmatrix}\begin{bmatrix} 2 & 1 \end{bmatrix}\\ & = \begin{bmatrix} \frac{4}{5}&\frac{2}{5} \\ \frac{2}{5}&\frac{1}{5} \end{bmatrix}\\ \Rightarrow \mathbf y = P\mathbf b & = \begin{bmatrix} \frac{4}{5}&\frac{2}{5} \\ \frac{2}{5}&\frac{1}{5} \end{bmatrix}\begin{bmatrix} -1 \\ 2 \end{bmatrix}\\ & = \begin{bmatrix} \frac{-4}{5}+\frac{4}{5} \\ \frac{-2}{5}+\frac{2}{5} \end{bmatrix}\\ \therefore P\mathbf b & = \begin{bmatrix} 0\\0 \end{bmatrix} = \mathbf 0\\ \end{align*}$$

To find $\delta \mathbf b$, will let $\delta \mathbf b = \mathbf b + \mathbf u$ where $\mathbf b = \begin{bmatrix} -1 \\ 2 \end{bmatrix}$ and $\mathbf u = \begin{bmatrix} x \\ y \end{bmatrix}$such that $P\delta \mathbf b = \delta \mathbf b$ that is; $$\begin{align*} \Rightarrow P(\mathbf b + \mathbf u) &= (\mathbf b + \mathbf u) \\ \Rightarrow \frac{1}{5} \begin{bmatrix} 4&2 \\ 2&1\end{bmatrix} \begin{bmatrix} -1+x \\ 2+ y \end{bmatrix} &= \bigg(\begin{bmatrix} -1+x \\ 2+ y \end{bmatrix}\bigg) \\\\ \Rightarrow \begin{bmatrix} 4x+2y \\ 2x+y &\end{bmatrix} &= \bigg(\begin{bmatrix} -5+5x \\ 10+ 5y \end{bmatrix}\bigg) \\ \end{align*}$$ Equating elements in the two vectors we obtain; $$\begin{align} x = 2y+5\\ -4y = 10-2x\end{align}$$ Subbstituting for x in equation (26), we obtain; $$\begin{align} -4y = 10-2(2y+5)\\ \Rightarrow 0= 0\end{align}$$ This means that $y$ can take on any real number ie. $y = t\in\mathbb R$.

Chosing $y = t =0$, this gives $x = 5$ and $\mathbf u = \begin{bmatrix} 5 \\ 0 \end{bmatrix}$

$$\begin{align*} \Rightarrow \delta \mathbf b &= \begin{bmatrix} -1+5 \\ 2+ 0 \end{bmatrix}\\ &= \begin{bmatrix} 4 \\ 2\end{bmatrix}\\ \therefore \delta \mathbf y = P \delta \mathbf b &= \frac{1}{5} \begin{bmatrix} 4&2 \\ 2&1\end{bmatrix}\begin{bmatrix} 4 \\ 2\end{bmatrix} = \begin{bmatrix} 4 \\ 2\end{bmatrix} = \delta \mathbf b \end{align*}$$

This can be geometrically illustrated as shown in the image below;

From the graph above, we see that the projection of $\mathbf b$ onto the $range(A)$ is a zero vector, that is $\mathbf y = P\mathbf b = \mathbf 0$. Surprisingly, a small change in $\mathbf b$ say $\delta \mathbf b$ results into non-zero pertubation $\delta \mathbf y$ in the projection. This means that the residual $\mathbf r = \mathbf b - A \mathbf x = \mathbf b$ and that the change in solution is very very large.

With reference to the above ideas, we conclude that when $\theta = \pi/2$, the condition number blows up i.e. $\kappa = \infty$.

The same conclusion can be reached by using the definition of the condition number of $\mathbf y$ with respect to perburbations in $\mathbf b$; $$\kappa=\frac{1}{cos(\theta)} = \frac{1}{cos(\pi/2)} = \infty$$

Problem 1(c)¶

Now consider the problem

$$\begin{equation} \begin{bmatrix} 2 \\ 1 \end{bmatrix} \begin{bmatrix} x \end{bmatrix} = \begin{bmatrix} 2 \\ 1 \end{bmatrix} \end{equation}$$

For this problem, show that $\kappa = 1$. What is qualitatively different about this problem than the problem in which $\kappa = \infty$?

Solution¶

Since in 1(a) we proved that $\kappa = \frac{1}{\cos(\theta)}$ and the angle between $\mathbf b$ and the $range(A)$ is zero, the condition numder for this case is then given as; $$\kappa=\frac{1}{cos(\theta)} = \frac{1}{cos(0)} = 1$$

For this problem, we see that the projection of $\mathbf b$ onto the column space of $A$ is $\mathbf b$ itself, that is $\mathbf y = P\mathbf b = \mathbf b$. This means that $\mathbf b$ is in the $range(A)$, the residual $\mathbf r = \mathbf b - A \mathbf x = \mathbf 0$ and that the change in the solution is very very small which isn't the case when $\kappa = \infty$.

Problem 2¶

---

Problem 18.1 in TB (page 136)

2(a)¶

2(b)¶

2(c)¶

2(d)¶

2(e)¶

(i)¶

The first condition number $K_{by}$ is attained when $\Vert\delta \mathbf b \Vert = \Vert P \delta \mathbf b\Vert$ and this occurs when $\delta\mathbf b$ is equal to zero except in the first $n = 2$ entries ($n$ is the number of columns of $A$),

I therefore used $\delta\mathbf b$ as $\begin{bmatrix} 0.0001 \\ 0.0001 \\ 0 \end{bmatrix}$ and this gave $K_{by} = 1.29108$

(ii)¶

The condition number $K_{bx}$ is attained when $\Vert A^\dagger \delta \mathbf b \Vert _2 = \Vert A^\dagger \Vert _2\Vert \delta \mathbf b\Vert _2$ and this occurs when $\delta\mathbf b$ is equal to zero except in the $2^{nd}$ entry.

For this I used $\delta\mathbf b = \begin{bmatrix} 0 \\ 10^{-8} \\ 0 \end{bmatrix}$ this gave $K_{by} = 54775.17695228298$

(iii)¶

The condition number $K_{Ay}$ is attained when $\delta A = \delta p v_2^T$ where $\delta p$ is the vector orthogonal to the $range (A)$ and $v_2$ is the second column of the right singular matrix $v$.

for this I used $\delta p = (I-A A^\dagger)\begin{bmatrix} 10^{-7} \\ 10^{-7} \\ 0 \end{bmatrix}$, $v_2^T = \begin{bmatrix} -0.70713035 &-0.70708321 \end{bmatrix}$ and obtained $\delta A = 54775.16415665296$

Problem 3¶

---

Show that if $(\lambda,\mathbf v)$ is an eigenvalue/eigenvector pair for matrix $A$, then $((\lambda-\mu)^{-1}, \mathbf v)$ is an eigenvalue/eigenvector pair for the matrix $(A - \mu I)^{-1}$.

Why is this observation useful when using the power iteration to find an eigenvalue close to $\mu$?

Solution¶

If $\lambda$ is an eigenvalue of matrix $A$ with a corresponding eigenvector $\mathbf v$, then we have that;

$$\begin{equation} A \mathbf v = \lambda \mathbf v \end{equation}$$

Subtracting $\mu \mathbf v$ from equation() above we obtain;

$$\begin{align} A \mathbf v - \mu \mathbf v = \lambda \mathbf v - \mu \mathbf v\\ (A-\mu I)\mathbf v = (\lambda-\mu)\mathbf v\\ %\Rightarrow \frac{1}{} = \frac{}{} \end{align}$$

Since $\mu \in \mathbb{R}$ is not an eigenvalue of $A$, then $A-\mu I$ is invertible and therefore multiplying the above equation with $(A-\mu I)^{-1}$, we get;

$$\begin{align} (A-\mu I)^{-1}(A-\mu I)\mathbf v &= (\lambda-\mu)(A-\mu I)^{-1} \mathbf v\\\\ \Rightarrow \mathbf v &= (\lambda-\mu) (A-\mu I)^{-1}\mathbf v\\\\ \Rightarrow (A-\mu I)^{-1}\mathbf v &= (\lambda-\mu)^{-1} \mathbf v \end{align}$$

Therefore from the above equation, we can conclude that $((\lambda-\mu)^{-1}, \mathbf v)$ is an eigenvalue/eigenvector pair for the matrix $(A - \mu I)^{-1}$.

Why is this observation useful when using the power iteration to find an eigenvalue close to $\mu$?

This observation is useful when using the power iteration to find an eigenvalue $(\lambda _j)$ close to $\mu$ because the dominant eigenvalue of $(A - \mu I)^{-1}$ will be $(\lambda _j -\mu)^{-1}$ from which we can easily obtain $\lambda _j$.

Problem 4¶

---

Exercise 29.1 (Lecture 29, TB page 223). This is a five part problem that asks you to code an eigenvalue solver for a real, symmetric matrix using the shifted $QR$ algorithm. Do your code in Python, using the Numpy qr algorithm where needed.

The basic steps are :

1. Reduce your matrix $A$ to tridiagonal form. You may use the hessenberg code we wrote in class.

2. Implement the unshifted $QR$ code (also done in class). Use the Numpy routine qr. Your iteration should stop when the off diagonal elements are smaller (in absolute value) than $\tau \approx 10^{-12}$.

3. Find all eigenvalues of a matrix $A$ using the "deflation" idea described in Algorithm 28.2.

4. Introduce the Wilkinson shift, described in Lecture 29.

Notes¶

• Your code should work for a real, symmetrix matrix

• Your code does not have to be efficient in the sense of optimizing the cost of matrix/vector multiplies and so on.

• Apply your algorithm to the Hilbert matrix scipy.linalg.hilbert. The entries of the $m \times m$ Hilbert matrix are given by

$$\begin{equation} H_{ij} = \frac{1}{i + j - 1}, \qquad i,j = 1,2,\dots m \end{equation}$$

Part (e)¶